function h_hat = predict_channel(y, simPar, predict_Method, user)
    
    if predict_Method == "LS"
        h_hat = y;
    elseif predict_Method == "MMSE_ISO"
%         mu = [0; 0; 10];% mu = [0.0001; -0.0001; 0.0001];
        mu = [0.0001; 0; 0.0001];
        theta = struct('rho', 2, 'v', [0; 0; 0], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, theta, simPar.k0, simPar.c);
%         eta = simPar.lambda/2/pi;
% n_azi = ( -(simPar.N-1)/2:1:(simPar.N-1)/2 )';
% % n_ele = ( -(N2-1)/2:1:(N2-1)/2 )';
% kernel_azi = exp(-0.1*pi*abs(simPar.BS.spacing*n_azi - simPar.BS.spacing*n_azi').^2);
% kernel_azi = besselj(0,0.85*abs(simPar.BS.spacing*n_azi - simPar.BS.spacing*n_azi')/eta);
% kernel_ele = exp(-0.096*pi*abs(d*n_ele - d*n_ele').^2/eta^2);
% Kernel_exp = kron(kernel_azi, kernel_azi);

%         h_hat =  kernel_azi*((kernel_azi + simPar.sigma2_n*eye(simPar.N))\y(33:64));
        h_hat =  Khb*((Kha + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

    elseif predict_Method == "GPR"
        if simPar.channelEstimator.GPR.mixedKernel == false
            initialTheta = struct( 'v',    [0;0;0], ...
                                'mu',   [0.001; 0; 0.001]);%[-0.2;+0.000; -0.1]
        
            [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
           [optKh, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta, simPar.k0, simPar.c);
            [optKhb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, optTheta, simPar.k0, simPar.c);

        else
            initialTheta = struct( 'v',    [0;0;1], ...
                                'mu',   [0.001; 0; 0.001]);
            [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
            optKh = zeros(simPar.N*simPar.T);
            optKhb = zeros(simPar.N,simPar.N*simPar.T);
            for idx = 1:simPar.channelEstimator.GPR.mixKernelNum
                [K, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta.kernelParams(idx), simPar.k0, simPar.c);
                [Kb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, optTheta.kernelParams(idx), simPar.k0, simPar.c);
                optKh = optKh + optTheta.kernelWeights(idx) * K;
                optKhb = optKhb+optTheta.kernelWeights(idx) * Kb;
            end
        end
        h_hat = optKhb*((optKh + simPar.sigma2_n*eye(simPar.N*simPar.T))\y); 
    
    elseif predict_Method == "GPRmu"
        mu = 1*user;
        theta = struct('rho', 2, 'v', [0; 0; simPar.v], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, theta, simPar.k0, simPar.c);
        h_hat =  Khb*((Kha + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

%         if simPar.channelEstimator.GPR.mixedKernel == false
%             initialTheta = struct( 'v',    [0;0;20], ...
%                                 'mu',   [20; 0; 5]);%[-0.2;+0.000; -0.1]
%         
%             [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
%             [optKh, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,optTheta, simPar.k0, simPar.c);
%             [optKhb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, optTheta, simPar.k0, simPar.c);
% 
%         else
%             initialTheta = struct( 'v',    [0;0;1], ...
%                                 'mu',   user);
%             [~, optTheta] = fit_GPR(y, initialTheta, simPar, false);
%             optKh = zeros(simPar.N);
%             for idx = 1:simPar.channelEstimator.GPR.mixKernelNum
%                 [K, ~] = get_K(simPar.x,simPar.x, simPar.p, optTheta.kernelParams(idx), simPar.k0, simPar.c);
%                 optKh = optKh + optTheta.kernelWeights(idx) * K;
%             end
%         end
%         h_hat = optKhb*((optKh + simPar.sigma2_n*eye(simPar.N*simPar.T))\y);

    elseif predict_Method == "AMP"
        A = dftmtx(simPar.N)/sqrt(simPar.N); 
        
        x_hat = camp2(y/sqrt(simPar.N), A, simPar.channelEstimator.AMP.lambda);
        h_hat = A*x_hat*sqrt(simPar.N); 
    
    elseif predict_Method == "OMP"
        A = dftmtx(simPar.N); 
        
        x_hat = OMP(y, A, 7);
        h_hat = A*x_hat; 

    elseif predict_Method == "GPRbes"
        [h_hat, ~] = fit_GPRbes(y, simPar); 
    
    elseif predict_Method == "Prony"
        R = y(:,1:simPar.T-1)'*y(:,1:simPar.T-1)/(simPar.N);
        v = y(:,1:simPar.T-1)'*y(:,simPar.T)/(simPar.N);
        P2 = inv(R)*v;
        p = pinv(y(:,1:simPar.T-1))*y(:,simPar.T);
        h_hat = y(:,2:simPar.T)*p;

    elseif predict_Method == "AR"
        fmax = simPar.v/simPar.lambda;
        T = simPar.Ts;
        z = 2*pi*fmax*T;
%         %一阶
%         aa = besselj(0, z);
%         h_hat = aa*y;
%         %二阶
%         Rxx=[1,besselj(0, z);besselj(0, z),1];
%         v_ = [besselj(0, z);besselj(0, 2*z)];
%         a_hat = inv(Rxx)*v_;
%         h_hat = a_hat(1)*y(1:simPar.N)+a_hat(2)*y(simPar.N+1:2*simPar.N);
        %T阶
        Rxx = toeplitz(besselj(0, (0:simPar.T-1)*z));
        Rxx = toeplitz(exp(-(0:simPar.T-1)*z));
        mu = [0.0001; 0; 0.0001];
        theta = struct('rho', 2, 'v', [0; 0; 20], 'mu', mu);
        [Kha, ~] = get_K(simPar.x,simPar.x, simPar.p, simPar.p,theta, simPar.k0, simPar.c);
        [Khb, ~] = get_K(simPar.x,simPar.xb, simPar.p,simPar.pb, theta, simPar.k0, simPar.c);
%         Rxx = zeros(simPar.T,simPar.T,'double');
        
        h_hat = zeros(simPar.N, 1);
        for jd=1:simPar.N
        for id=1:simPar.T
            r = zeros(1,simPar.T,'double');
            v_ = zeros(1,simPar.T,'double');
            r(id) = Kha(jd,(id-1)*simPar.N+jd);
            v_(id) = Khb(jd,(id-1)*simPar.N+jd);
        end
%         Rxx = toeplitz(r);
        y2 = [y(1:simPar.N),y(simPar.N+1:simPar.N*2)];
%         Rxx = y(1:simPar.N)'*y(1:simPar.N);

        v_ = [besselj(0, (1:simPar.T)*z)];
        v_ = [exp(-(1:simPar.T)*z)];
        v_ = v_.';
%         v_ = y(1:simPar.N)'*y(simPar.N+1:simPar.N*2);

%         a_hat = Rxx\v_;
        a_hat = pinv(Rxx)*v_;
        
        
simPar.sta = 350;
    R = zeros(simPar.T,simPar.T);
    DD = zeros(simPar.T+1,simPar.sta);
    v = zeros(simPar.T,1);
    simPar.a = zeros(simPar.N,simPar.T);
    D = zeros(simPar.N,simPar.T+1,simPar.sta);
    Num = zeros(simPar.T+1,simPar.sta);
        cdl_model = nrCDLChannel('ChannelFiltering', false);
        cdl_model.MaximumDopplerShift       = 2*simPar.v/simPar.lambda;
        cdl_model.UTDirectionOfTravel       = [0;0];
        cdl_model.CarrierFrequency          = simPar.fc;
        cdl_model.SampleRate                = 1/simPar.Ts;
        cdl_model.SampleDensity             = Inf;
        cdl_model.NumTimeSamples            = simPar.T+simPar.sta;
        cdl_model.TransmitAntennaArray.Size = [simPar.N, 1, 1, 1, 1];   % [M N P Mg Ng]
        cdl_model.ReceiveAntennaArray.Size  = [1, 1, 1, 1, 1]; 
        cdl_model.Seed                      = randi(10000);
        [pathGains, ~] = cdl_model();
        for idx1 = 1:simPar.T+simPar.sta
             h1(:,idx1) = reshape(sum(pathGains(idx1, :, :, :), 2), [simPar.N, 1]);
        end
        for j = 1:simPar.T+1
        simPar.h((j-1)*simPar.N+1:j*simPar.N) = h1(simPar.sta+j-1);
        end
        for idx2 = 1:simPar.T+1
            for idx3 = 1:simPar.sta
           D(:,idx2,idx3) = h1(:,idx3+idx2-1);
           Num(idx2,idx3) = idx3+idx2-1;
            end
        end
    for idx = 1:simPar.N
    DD(:,:) = D(idx,:,:);
    RQ = DD*DD'/simPar.sta;
    R = RQ(1:simPar.T,1:simPar.T)+eye(simPar.T)*0.01;
    v = flip(RQ((simPar.T+1),(1:(simPar.T))));
    v = v.';
%     [V,D2,W] = eig(R)
    simPar.a(idx,:) = inv(R)*v;
    end
        for i = 1:simPar.T
%             h_hat = h_hat+a_hat(i)*y((simPar.T-i)*simPar.N+1:(simPar.T-i+1)*simPar.N);
%             h_hat(jd) = h_hat(jd)+a_hat(i)*y((simPar.T-i)*simPar.N+jd);
            h_hat(jd) = h_hat(jd)+simPar.a(jd,i)*y((simPar.T-i)*simPar.N+jd);
        end

        end
%  v_ = [exp(-(1:simPar.T)*z)];
%         v_ = v_.';
% 
% %         a_hat = Rxx\v_;
%         a_hat = inv(Rxx)*v_;
%         h_hat = zeros(simPar.N, 1);
% 
%         for i = 1:simPar.T
%             h_hat = h_hat+a_hat(i)*y((simPar.T-i)*simPar.N+1:(simPar.T-i+1)*simPar.N);
%         end
%         a_hat = zeros(simPar.N,2);
%         h_hat = zeros(simPar.N,1);
%         p = 2;%order
%         for i=1:simPar.N
%             Rxx = zeros(p, p); % 自相关矩阵
%             v = zeros(p, 1); % Yule-Walker方程的向量
%             % 计算自相关矩阵和Yule-Walker方程的向量
%             for j = 1:p
%                 for k = 1:p
%                 if k == j
%                     Rxx(lag+1, lag+1) = abs(y(i)).^2;
%                 else
%                     m = 0;
%                     for ii = 1:p
%                         m = m+abs(y(i+(k-1)*simPar.N)).*conj(H(:, lag+1:end));
%                     end   
%                     Rxx(j, k) = sum(abs(y(i+(k-1)*simPar.N)).*conj(H(:, lag+1:end)))/p;
%                 end
%                 end
%             end
%                 for k = 1:p
%                     if lag == 0
%                         v(k) = sum(abs(H(:, k)).*conj(H(:, 1)));
%                     else
%                         v(k) = sum(sum(abs(H(:, k+lag:end)).*conj(H(:, 1:end-lag))));
%                     end
%                 end
%             
%         end

%         for i=1:simPar.N
%             [a_hat(i,:)] = ar([y(i),y(simPar.N+i)],2,'ls');
%             h_hat(i) = a_hat(i,1)*y(i)+a_hat(i,2)*y(simPar.N+i);
%         end

    else
        error("Error: Channel estimation method %s is not implemented.", predict_Method);
    end

end

